{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load the required packages\n",
    "using ODE\n",
    "using JLD\n",
    "using ForwardDiff\n",
    "\n",
    "set_bigfloat_precision(113)\n",
    "\n",
    "# Define the system for the solver\n",
    "function roberAD(x)\n",
    "    k1 = parse(BigFloat,\"0.04\");\n",
    "    k2 = parse(BigFloat,\"1e4\");\n",
    "    k3 = parse(BigFloat,\"3e7\");\n",
    "    \n",
    "    return [-k1*x[1]+k2*x[2]*x[3],\n",
    "    k1*x[1]-k2*x[2]*x[3]-k3*(x[2])^2,\n",
    "    k3*(x[2])^2]\n",
    "end\n",
    "\n",
    "function rober(t,x)\n",
    "    return roberAD(x)\n",
    "end\n",
    "\n",
    "function getJacobian(t,x)\n",
    "    return ForwardDiff.jacobian(roberAD,x);\n",
    "end\n",
    "\n",
    "# Set up the initial conditions\n",
    "tSpan = [zero(BigFloat);parse(BigFloat,\"10.0\").^collect(0:11)];\n",
    "x0 = [one(BigFloat),zero(BigFloat),zero(BigFloat)];\n",
    "\n",
    "# Set the tolerances\n",
    "# ATol = RTol*1e-6\n",
    "RTol = parse(BigFloat,\"1e-20\");\n",
    "ATol = parse(BigFloat,\"1e-26\");\n",
    "\n",
    "# Solve and get the solution at T = tEnd\n",
    "(t,x_tmp) = ode23s(rober,x0,tSpan;\n",
    "reltol=RTol,abstol=ATol,points=:specified,\n",
    "jacobian = getJacobian,\n",
    "minstep=parse(BigFloat,\"1e-8\"));\n",
    "\n",
    "x_ref = Array{BigFloat}(11,3);\n",
    "\n",
    "for i=1:11\n",
    "    x_ref[i,:] = x_tmp[i+1,1][:];\n",
    "end\n",
    "\n",
    "# Save the solution to a file\n",
    "save(\"refSolRober.jld\",\"x_ref\",x_ref);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
